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ABSTRACT 

We present the results of a photometric monitoring campaign of three well stud- 
ied FU Orionis systems (FU Orionis, V1057 Cygni and V1515 Cygni) undertaken at 
Maidanak Observatory between 1981 and 2003. When combined with photometric 
data in the literature, this database provides a valuable resource for searching for 
short timescale variability - both periodic and aperiodic - as well as for studying the 
secular evolution of these systems. In the case of V1057 Cyg (which is the system 
exhibiting the largest changes in brightness since it went into outburst) we compare 
the photometric data with time dependent models. We show that prior to the end of 
the 'plateau' stage in 1996, the evolution of V1057 Cyg in the V — (B — V) colour- 
magnitude diagram is well represented by disc instability models in which the outburst 
is triggered by some agent - such as an orbiting planet - in the inner disc. Following the 
end of the plateau phase in 1996, the dimming and irregular variations are consistent 
with occultation of the source by a variable dust screen, which has previously been 
interpreted in terms of dust condensation events in the observed disc wind. Here we 
instead suggest that this effect results from the interaction between the wind and an 
infalling dusty envelope, the existence of this envelope having been previously invoked 
in order to explain the mid infrared emission of FU Orionis systems. We discuss how 
this model may explain some of the photometric and spectroscopic characteristics of 
FU Orionis systems in general. 
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1 INTRODUCTION 

Despite being a rather small class of young stellar objects, 
the outbursting FU Orionis systems have been studied ex- 
tensively over the last 20 years, since they represent an 
interesting "laboratory" to study the accretion process on 
newly born stars and the interaction between the protostar 
and its environmen t. In fact, according to the most pop- 
ular interpretation jHartmann fc KenvonllT'996ll . their out- 
bursts are associated with episodes of enhanced accretion 
through an other wise "normal", T T auri like, protostellar 
disc (see, however. iHerbig et a l. 2003 for a dissenting view). 
Support to this interpretation comes mainly from the mod- 
eling of the broad band Spectral Energy Distribution (SED) 
in terms of accretion disc SEDs, from their apparent spectral 
type, which is characterized by high effective temperatures 
at smaller wavelength and by lower effective temperatures at 
larger wavelengths, and from the shape of some optical and 
near-infrared absorption lines, which show the typical double 



peak expected from a rotating disc lIHartmann fc Kenvonl 
ll985lll987l:lKenvon et allll98ai . 

The time evolution of the outbursts may vary consid- 
erably from object to object. If we refer to the three best 
studied objects (FU Ori itself, V1057 Cyg and V1515 Cyg, 
the only ones for which a detailed light curve has been 
recorded since the beginning of the outbursts), it can be 
seen that while for some of them (FU Ori and V1057 Cyg) 
the transition to the outburst phase is rather abrupt (with 
rise timescale of the order of 1 yr), for V1515 Cyg it is con- 
siderably more gentle (rise timescale m 20 yrs). The post- 
outburst light curve is also quite different for different ob- 
jects: while FU Ori has remained almost steadily in the 
high luminosity state for at least 70 yrs, V1057 Cyg has 
faded rapidly (with a timescale of roughly 10 yrs) and then 
showed a "plateau" in the light curve, until eventually un- 
dergoing an abrupt luminosity drop in the mid 90s, followed 
by a rather erratic photometric variability. V1515 Cyg, on 
the other hand (which is characterized by a smaller peak 
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luminosity with respect to the other two), while remaining 
overall in the high luminosity state, has often shown strong 
photometric variability. It would then appear to be difficult 
to reconcile this wide variety of behaviours within a single 
model. 

However, the differences in rise timescales can be un- 
derstood when one considers the detailed outburst mecha- 
nism. This is generally considered to be due to the onset 
of a thermal instability in a p rotostellar disc fed at a high 
enough rate jBell fc Linl 11994) . a mechanism analogous to 
the one considered in dwarf novae outbursts. In this pic- 
ture, partial ionization of hydrogen in the inner disc gives 
rise to a thermal instability and to a limit cycle behaviour, 
that ca n be able to rep roduce repetitive outbursts. It can be 
shown <Lin et alJll985l) that outbursts initiated at the inner 
edge of the disc propagate inside-out, resulting in a long rise 
timescale, while outbursts triggered somehow far from the 
inner edge are characterized by short rise timescales. It is 
then possible to associate the long rise timescale of V1515 
Cyg to a non-triggered outbursts, and the fast rise in V1057 
and FU Ori to a triggered one. Models of triggered outbursts 
have been constructed, both in the case where the outburs t 
is artificially triggered ijClarke et alJll99ft IrBell et al]ll99it . 
and where the trigger is provided by the inter action between 
the disc and an embedded massive planet dClarke fc Sverl 
ll996l:lLodato fc Clarkdliooi) . 

The environment of FU Orionis objects shows further 
complexities: there is, in fact, substantial evidence for the 
presence of an infalling envelope, a bipolar jet and strong 
mass loss in the form of a disc wind. 

The best evidence for the presence of an envelope is rep- 
resented by the mid and far infrared SED, that shows a sig- 
nificant excess with r espect to the predictions of s tandard ac- 
cretion disc models jKenvon fc Hartmannlll99ll henceforth 
KH91). In many cases, even the inclusion of reprocessing of 
the inner disc luminosity by an outer, flared disc is not suf- 
ficient, because it would require too large degree of flaring. 
Even though different possible explanati ons for this long- 
wave l ength excess have be en discussed llLodato fc Bertinl 
l200ll : lAbraham et alJ l2004f) , it is commonly attributed to 
reprocessing from a dusty infalling envelope, whose inner 
radius should be located at ~ 10 au from the central source 
JKH91I) . This model however, requires the presence of a cav- 
ity in the envelope through which the central optical source 
is seen. This cavity can be produced by a bipolar jet (typi- 
cally observed also in many T Tauri stars). The presence of 
spatially separated knots in the jets of many Herbig-Haro 
objects has been indeed interpreted as d ue to recurrent FU 
Orion is activity in young stellar objects llReipurth fc Aspml 
I1997D . 

Strong winds have been observed in FU Orionis objects 
using a variety of different tracers: optical and near-infrared 
lines dCroswell et al.lll987l) . that show the characteristic P 
Cygni profi le, lower excitation lines o f neutral metals and 
TiO bands iHartma m]_j£j<gny^3n|l99 ^) and even c ontinuum 
radio observations lRwjTrgu^z^^HartmandlT992l) . Typical 
wind velocities are observed to be in the range 300 — 400 
km/sec and mass loss rates are observed to vary between dif- 
ferent objects, being very large for FU Ori (« 10 -5 MQ/yr) 
and roughly one order of magnitude smaller for V1057 Cyg 
during the plateau phase and V1515 Cyg. 

In addition, in some cases there is also evidence for 



the presence of comp anions: the FU Orio nis object Z CMa 
is a binary system llKoresko et alJ Il991 ), and recently a 
companion to FU Ori has been found (IWang et al ] 12004 
iReipurth fc AspinlEooi) . Note, however, that this compan- 
io n is unlikely to be the cau se of the outburst, as suggested 
bv lBonnell fc BastienUl992D . since its separation (f» 200 au) 
is too large to reproduce the fast rise of the outburst. 

In this paper we present the result of a photometric 
monitoring campaign of three well studied FU Orionis sys- 
tems (FU Orionis, V1057 Cygni and V1515 Cygni) under- 
taken at Maidanak Observatory between 1981 and 2003. 
These new data, combined with historical photometric data, 
are then used to provide a quantitative test to theoretical 
models. 

In particular, we use a thermal instability model to de- 
scribe the onset of the outburst and compare the colour evo- 
lution predicted by this model to the observations of V1057 
Cyg and V1515 Cyg. We obtain a good match to the ob- 
servations, better than what can be obtained using a series 
of stead y state accretion disc models, as done in the past 
JKH91I) . As mentioned above, the long rise time and the 
colour evolution of V1515 Cyg are well described by non- 
triggered thermal instability models, whereas the data for 
V1057 Cyg require a triggered model, where the triggering 
mechanism is here taken to be the interaction with a small 
mass companion in the disc llLodato fc Clarkell2004D . 

Furthermore, we present a new model of the interaction 
between the disc wind and the infalling dusty envelope. The 
main assumption of this model is that the wind strength 
scales with the luminosity of the disc. We find that suffi- 
ciently strong winds (like that of FU Ori) are able to blow 
the envelope out to large distances, while weaker winds (like 
the one of V1515 Cyg) do not. For V1057 Cyg, the wind 
initially is able to blow the envelope away to large distance, 
but as the luminosity of the disc decreases the envelope falls 
back to small radii over a timescale of « 20 yrs. The dusty 
envelope provides occultation of the central source. In this 
way, we are then able to explain the diverse post-outburst 
photometric variability of the three objects within a single 
simple model. 

The paper is organised as follows. In section [5] we de- 
scribe our new observations. In section |3] we show the com- 
parison between the observed colour evolution of V1057 Cyg 
and V1515 Cyg and the one predicted by disc thermal in- 
stability models. In section 2] we describe our model for the 
wind-envelope interaction and in section|K|we draw our con- 
clusions. 



2 OBSERVATIONS 
2.1 New photometry 

The new UBVR photometric data for FU Ori, V1057 Cyg 
and V1515 Cyg presented here were taken at Maidanak Ob- 
servatory during a long monitoring campaign from 1981 (for 
FU Ori, from 1984) to 2003, using two 60-cm Zeiss reflec- 
tors and the 48-cm AZT-14 reflector equipped with iden- 
tical photon-counting photometers, which reproduced the 
Morgan- Johnson system. The observations were carried out 
by differential method using comparison stars. The rms er- 
ror of a single measurement was AV — A(V —R) = m .015, 
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Figure 1. Light-curves of V1057 Cyg (left ), V1515 Cyg (middle) and FU Ori (right). Blue squares ind ic ate our new data, wh i le red 
trian g les indicate historical data taken from | Mendoza| il97%lRieke et alj<1972l): ISchwartz fc Snowl il972fr : lBossei] il97^ : rWelinl <197Et 



ll976ft : lLandoltj h975lll977ft ; lKolotilovl il977ft : lKopatskavjll98^ : llbragimov fc Shevchenkd il98- 



A(B - V) = m .02. Only for FU Ori itself the rms error 
in U — B was about m .03, whereas for V1057 Cyg and 
for V1515 Cyg this one was much higher because of their 
brightness in U and we used measurements at U — B only 
for approximate estimates. 

Full tables with the detailed photometry for the three 
program objects are available electronica l ly. Part of these 
data were already published in llbrahimovl il996tlT999l) . 



2.2 Long-term light curves 

Fig- shows the long-term photometry of the three objects 
considered combining historical data (red triangles, see cap- 
tion for references) and our new data (blue squares). 

Among the three objects considered, the most interest- 
ing long-term behaviour is the one of V1057 Cyg. This star, 
in fact, is the one which shows the fastest time evolution: 
from the peak of the outburst in the early '70s to 2003 its 
luminosity has decreased by almost 5 magnitudes in B and 
by 4 magnitudes in V . Its photometric behaviour after the 
outburst can be divided in three phases: a period of rather 
steep decrease of brightness, a "plateau" phase, and a pe- 
riod of smooth variability after a sharp drop in luminosity 
in 1995. Maidanak photometry provides a good coverage of 
these two last parts of the light curve. The long-term strong 
evolution of the system makes it an ideal candidate to test 
evolutionary outburst models. We have performed this test 
by comparing theoretical models to the colour-magnitude 
V — (B — V) diagrams (Fig. |7J , discussed below in Section 

m 

After the luminosity drop in 1995, slow year-to-year 
variability was observed in V1057 (see Fig. 0. Apart from 
this variability, V1057 also shows smooth variations within 
one year, but while during the "plateau" phase the max- 
imum amplitude of this variability was about m .2, it in- 
creased to m .5 after the luminosity drop in 1995. We sub- 
tracted long-term variability from the lightcurve between 
1995 and 2003 and analysed the residual time series us- 
ing the Starlink CL EAN algorithm (originally developed by 
iRoberts et allll987l) to look for periodicity in the short-term 



variations. Although we detected a period of 150 days dur- 
ing 1995-1999, this period was not present in each of the 
individual seasons and is, in fact, uncomfortably close to 
the total duration of each of the seasonal light curves. 

V1057 Cyg shows also short-term variability and we 
analysed the time series to look for short-term periodicity (of 
order several days) and to test for flickering behaviour. We 
discuss in more detail this variability in Section [2.31 below . 

V1515 Cyg did not show any sharp decrease in luminos- 
ity during the last years, but we can see a gradual decrease 
in brightness during 1996-1999 (see Fig. 0. Short-term lu- 
min osity drops have be en observed in 1980 and 1987 (see 
also iKenvon etai]|l99lf) . 

During 2000-2003 FU Ori continued a slow decreasing 
of its brightness (Fig.0. FU Ori also did not show any sharp 
luminosity drops in recent years. 

For these last two objects, given their small luminosity 
evolution, a time-dependent theoretical model would be less 
meaningful than for V1057 Cyg. However, as an example, 
we have constructed a disc thermal instability model and 
compared it with observations also for V1515 Cyg and we 
describe it below in Section |3] 

2.3 Short-term variability and flickering 

IKenvon et alJ <l2000l) have obtained some evidence for flicker- 
ing - a non periodic brightness variability with an amplitude 
of roughly m .2, on timescales of the order of a few days - 
in FU Ori. We have analysed our data for V1057 Cyg and 
for V1515 Cyg to look for a similar flickering behaviour. 

2.3.1 VI 057 Cyg 

We have subtracted any long-term variation from each sea- 
sonal data set from 1986 to 1990 (i.e. during the photometric 
"plateau") by fitting a polynomial to the observed light- 
curve. We have then analysed all datasets with CLEAN. 
Only in 1989 (Fig. our analysis indicated the presence 
of a periodic component with high probability in residual 
light-curve for BVR data (false alarm probability smaller 
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Figure 2. BVR phased light-curves for V1057 Cyg in 1989. The 
period suggested by these plots in P « 14 days. 





Figure 3. Correlation between V magnitude and B — V (left) and 

V — R colours for V1057 Cyg in 1988, in which year no periodic- 
ity of the variability was revealed by our analysis. Brightness and 

V — R colour variations appear to correlate well. B — V colours 
appear to display a larger scatter, but some correlation can still 
be noted. Dotted lines show the direction of the interstellar red- 
dening vector. 
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Figure 4. Correlation between brightness and colour variations 
for V1515 Cyg in 1996, during which year no periodicity in the 
variability was revealed by our analysis. As for V1057 Cyg, V — R 
colour variations correlate well with brightness variations. Dotted 
lines show the direction of the interstellar reddening vector. 



than 1%). Periodograms suggest a period of 14 ± 0.2 days 
in all three bands, where the error represents the half width 
half maximum of the periodogram peak. 

For all other years, although we observe a variabil- 
ity with amplitude ~ m .l in V, no periodicity was de- 
tected, as expected for a flickering variability. To verify that 
these changes in bri ghtness and colours are real, we follow 
iKenvon et alJ i2000f) and plot B — V and V — R colours 
against V for 1988 (for which year no periodicity was de- 
tected). The variations on short timescales clearly exceed 
photometric errors and colour changes of V — R clearly cor- 
relate with brightness changes (Fig.|3J. Although the scatter 
in the V — (B — V) diagram is bigger than in V — (V — R), 
we can also see some c orrelation between V and B — V (note 
a similar result also in IKenvon et alll200Cl for FU Ori). Note 
also that the amplitude of variation in V 1 057 C yg is roughly 
a half of that detected by IKenvon et alJ fcOOCh for FU Ori. 



2.3.2 V1515 Cyg 

We performed a similar analysis also for V1515 Cyg. After 
subtracting any long-term evolution we observe a short-term 
variability with an amplitude ~ m .3 in V . Analysis with 
CLEAN only revealed some periodicity on short time-scales 
from the 1987 season (although with a slightly higher false 
alarm probability of 1.5 %). As for V1057 Cyg, non-periodic 
colour variations in V — R correlate well with brightness 
variations (Fig. • Also in this case the scatter in B — V 
is much larger, but we can still see some correlation in the 
V — (B — V) diagram. 

As mentioned above, our analysis indicates a periodic 
component with high enough probability in residual light 
curves for BVR data only for the 1987 season. Periodograms 
suggest a period of 13.9 ±0.1 days in all three bands (see 
phased light curves, Fig. 
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Figure 5. BVR phased light-curves for V1515 Cyg in 1987, the 
only year in which our analysis revealed some periodicity in the 
photometric variability. The period suggested is P = 13.89 days. 



2.4 Correlation between spectral and photometric 
variability 

FU Orionis objects show also a wealth of spectral variations 
on short timescales (of the or der of days) . These ar e reviewed 
and discussed extensively in IHerbig et alJ (I2003I) . They in- 
clude spectral variations of the equivalent width of the Ha 
profiles. For examp le, for FU Ori itself, during 1997-1999 
IHerbig et alJ i2003l) detected a variability with a period of 
roughly 14 days. This kind of variability h as als o been ob- 
served earlier in FU Ori by lErrico et alJ (J2003) , but they 
detected a smaller per iod of ~ 6 day s. 

Most noticeably, IHerbig et alJ (|2003Fl detect a peri- 
odic modulation of the double-peaked line profiles in the 
cross-correlation functions of optical photospheric absorp- 
tion lines. The periods detected in 1997 are s=s 3 days for 
FU Ori and ~ 4 days for V1057 Cyg (even if the latter 



Figure 6. Plot of our V magnitude p hotometry against s imulta- 
neous radial velocity measurements bv lHerbig et al.l J20031) , as de- 
rived from variability in the shape of double-peaked photospheric 
absorption lines. No apparent correlation is found. 



period is more uncertain') . IClarke fc Armitagd J2003) have 
shown how a periodic modulation of these lines can be due 
to the presence of a massive planet embedded in the FU Ori- 
onis discs. The observed periods in V1057 Cyg and FU Ori 
would then correspond to an orbit with a semi-major axis 
of roughly IORq, assuming that the planet is in Keplerian 
rotation around the host star. This is indeed the distance 
at which we expect to find a massive planet if it is respon- 
sible for triggering the thermal instability in the disc (see 
lLodato fc Clarkdl2003 and Section ^ below). Note that the 
two objects for which this modulation has been observed are 
fast-rise FU Orionis objects, for which a triggered outburst 
is needed. 

In the model of lClarke fc Armitaed i2003l) . however, the 
spectral variability should not be correlated with any pho- 
tometric variability. This is because the emission pattern is 
in this case constant in a frame co-rotating with the pu- 
tative planet and the modulation in the line profiles arises 
only from the changing location of the planet's orbital veloc- 
ity vector to the line of sight. We have then tested whether 
there is any correlation between our observed photometric 
va riability and the mo dulation of the line profiles observed 
bv lHerbig et al.l i2003t) for those periods where our obser- 
vations overlap with theirs. Fig. |S| shows the results of this 
analysis for V1057 Cyg, where our V luminos ity is plotted 
agains t coeval radial velocity measurements bv lHerbig et alJ 
(2003). No clear correlation is indeed found. A similar result 
holds also for the B band. 

Our observatio ns of FU Ori do not overlap with those of 
IHerbig et alJ J2003), but we used our photometry obtained 
from 1997 to 1999 to search for periods in the range of a few 
days on the light curves (i.e. with the same periods as the 
observed spectroscopic period). The analysis did not reveal 
any reliable photometric period in this epoch. 
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3 PHOTOMETRIC EVOLUTION: DISC 
THERMAL INSTABILITY MODELS 

In this section we describe the disc thermal instability mod- 
els that we have constructed in order to describe the colour 
evolution of the FU Orionis objects considered. We have con- 
sidered a triggered outburst for the fast-rise object V1057 
Cyg and an untriggered one for V1515 Cyg. We have not 
considered the case of FU Ori itself for two reasons: firstly its 
luminosity has changed very little after the outburst (total 
variation in B is approxima tely one magnitude), sec ondly, 
as discussed in more detail in lLodato fc Clarke! i2004) . ther- 
mal instability models are generally not able to reproduce 
such a long lived (more than 70 yrs) outburst. 



3.1 The colour evolution of V1057 Cyg 

V1057 Cyg is, among all the FU Orionis objects, the one 
that displays the fastest decay from the outburst phase. In 
fact, it has faded from the peak bolometric luminosity of 
8001/0 to « 2001/0 over a timescale of about 20 yrs. 

The time evolution of the Spectral Energy Distribution 
(SED) of V1057 C yg has been modeled in term s of an accre- 
tion disc SED bv lKenvon fc Hartmannl il99ll) . They com- 
pared the observed SED of V1057 Cyg at different epochs 
with synthetic disc models including reprocessing from an 
infalling envelope, in order to account for the observed far 
infrared excess. For each epoch, they assumed the disc to be 
in a steady state, characterized by a constant mass accretion 
rate M, whose value was chosen so as to match the instanta- 
neous luminosity of the system. They assumed a reddening 
correction of Av = 3.5 and found an overall good agree- 
ment with the data, except that the models tended to be 
significantly redder than the observations at early epochs, 
i.e. close to the luminosity maximum. However, these mod- 
els are clearly limited by the assumption that at any given 
epoch the disc is assumed to be in a steady state, whereas 
it is in fact evolving quite rapidly. In the following we will 
describe how a better agreement with observations of the 
colour evolution of V1057 Cyg can be achieved when prop- 
erly using time-dependent accretion disc models. 

We have constructed time dependent ID accretion disc 
models to account for the onset and evolution of the out- 
burst, in the framework of the thermal i nstability model 
JClarke et alJll989l 1199(3 : iBell fc Linlll994l) . Full details on 
the thermal instability model and on our numerical model- 
ing of the outburst can be found in lLodato fc Clarke! J2004) . 

VI057 Cyg is one of the fast rise FU Orionis objects. In 
fact, it shows a rise timescale of the order of one year, com- 
parable to the one of FU Ori, but much sh orter than the one 
of V1 515 Cyg (« 20 yrs). We have shown jLodato fc Clarke! 
120041) that a fast rising light curve can be obtained if the 
thermal instability, r ather than being initiated at the in- 
ner disc radius (as in lBell fc Linlll994 models), is triggered 
within the disc, at a radius of ~ 107?©. Triggered thermal 
instabilities also result in a much larger amplitude of the out- 
burst, with peak luminosity reaching 10 3 Lq (to be compared 
with peak lumi nosities of ~ 10 2 L© obt ained in non-triggered 
outbursts, see lLodato fc Clarke! 12004) . This compares well 
with the observations, since the fast rise FU Orionis objects 
(FU Ori and V1057 Cyg) have indeed a larger peak luminos- 
ity with respect to V1515 Cyg. It then appears reasonable 



to describe FU Ori and V1057 Cyg in terms of triggered 
thermal instability models, and V1515 Cyg in terms of non- 
triggered models. 

Rather than a ssuming an ad hoc triggering mecha nism, 
as previously done dClarke et al ll989UBell et al Jl995h . here 
we will considered the case where the instability is trig- 
gered by the presence of a small mass companion (a massive 
planet, for example) embedded in the disc. 

The details of the planet-disc interaction assumed here, 
and its inclusion in thermal in stability models can be 
found in lLodato fc Clarke! (l2004f) . Here we briefly summa- 
rize the basic physical mechanism responsible for the trig- 
gering. During quiescence, a massive planet (of the order of 
10Mj U piter) is able to open up a deep gap in the disc and will 
undergo Type II migration. However, in the case in which 
the mass of the planet is large compared to the local disc 
mass, 

M p > AtyE(R p )R 2 p (1) 

(where R p is the position of the planet and S(i? p ) is the 
local disc density), the actual migration timescale is longer 
than in normal Type II migration, and as a consequence disc 
material will bank up upstream of the planet, increasing the 
surface density of the disc and eventually leading to a ther- 
mal instability. We have already shown llLodato fc Clarke! 
2004) the effectiveness of this model in producing fast rising 
light curves in FU Orionis objects. 

The basic physical parameters that determine the be- 
haviour of the outburst are the mass of the planet M p and 
the input mass accretion rate at which the disc is fed by 
the envelope at large radii, M env . Additional parameters en- 
tering the models are the value of the viscosity parameter 
a, the mass of the central star M*, and the inner radius of 
the disc 7? m in- In the following we will take M 4 = 0.5Mq, 
Rmin = 5Rq and (as customary in thermal instability mod- 
els) we will assume two different values for a on the lower 
and upper branch of the thermally stable equilibrium curves: 
we have taken qi ow = 10 -4 and Qhigh = 10 -3 . We have 
then varied the value of M cnv and of M p in order to match 
the observed colour evolution of V1057 Cyg. The best fit 
models have M p = 10Mj U pit er and M cn v = 10 Mq/jt. 
For this choice of parameters the outburst is triggered at 
_R tr i g w 12Rq, and we obtain a fast-rise outburst with 
rise timescale of the order of ~ 2 yrs, with a peak lumi- 
nosity ipcak ~ 750-Lq and a peak mass accretion rate of 
Mpoak ~ 6 x 10~ 4 M Q /yr. 

At any given time during the outburst, the model pro- 
vides us with the full surface temperature profile of the disc, 
from which we can compute the broad band fluxes by sim- 
ply assuming a blackbody emission from every annulus of 
the disc. Our models also include the contribution of the 
central star. This contribution is, however, negligible when 
compared to the high luminosity arising from the disc. The 
disc is assumed to be viewed face-o n and the distance to 
V105 7 Cyg is assumed to be 600 pc llKenvon fc Hartmannl 
Il99ll) . Note that (here as in the case of V1515 Cyg, dis- 
cussed below) a small inclination of « 30°, required to fit 
the broadening of optical absorption lines, can be easily ac- 
comodated by changing the assu med distance within the un- 
certainties llKenvon et alj[T988l) . The fluxes have then been 
reddened assuming a redden ing correction of Av ~ 3.3 (con- 
sistent with the estimate of lBell et alj [l995). The resulting 
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Figure 7. Colour evolution of V1057 Cyg. The red triangles show 
the historical data points, while the blue squares show the new 
observations described in this work. The solid line shows the evo- 
lution of a thermal instability model triggered by a 10Mj up ; ter 
planet. The dashed line show the colour evolution predicted by a 
series of steady state models with different M. The short-dashed 
line indicates the direction of the extinction vector, for a standard 
interstellar extinction law. 



colour evolution in the V — (B — V) plane is shown in Fig. 
[7| with a solid line. The red triangles show the historical 
data, while the blue squares show our new observations, de- 
scribed above. The dashed line shows the colours of a se- 
ries of steady state models, with different mass accretion 
rates. These latter models a re similar to those described by 
iKenvon fc Hartmannl <199ll) . 

Note that, for a given V magnitude, close to the peak 
of the outburst the time dependent models are bluer than 
the corresponding steady state model. This is because in the 
time- dependent case, at the beginning of the outburst, while 
the unstable front propagates outwards, only the innermost 
part of the disc are in the high state, while the outer parts, 
which contribute to longer wavelengths, are still in the low 
state and therefore have a much lower luminosity with re- 
spect to steady state model, that assumes that all the disc 
is in outburst. In this way, our time dependent model is 
able to obtain a much better fit to the observations at all 
epochs during the outburst, while the series of steady states 
is redder than ob served close to the peak of th e outburst, as 
already noted bv lKenvon fc Hartm ann (1991). 

Another thing to notice is that the simple viscous disc 
evolution is not able to account for the recent drop in lu- 
minosity occurred in 1995. This is for two reasons: (i) the 
drop is too abrupt to be due to viscous evolution (see light 
curve in Fig.0, (ii) the colour change during the drop can- 
not be reproduced by the time-dependent models. In fact, 
all the data points after the drop occupy a different region 
in the V — (B — V) plane with respect to earlier observations 
(see Fig. Q . The colour variation observed during the lumi- 
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Figure 8. Colour evolution of V1515 Cyg. The red triangles show 
the historical data points, while the red squares show the new 
observations described here. The solid line shows the evolution of 
an untriggered thermal instability model. The dashed line show 
the colour evolution predicted by a series of steady state models 
with different M. 

nosity drop is much redder than the one predicted by the 
viscous evolution model. On the other hand the direction of 
the colour variation is consistent with a standard interstellar 
extinction law, adopting R — 3.3 (the direction of the extinc- 
tion vector is indicated with a short-dashed line in Fig. 0. 
We are therefore led to the conclusion that the recent drop in 
luminosity is due to a sudden enhancement of the extinction 
along our line of sight to VI 057 Cyg, rather than to the evo - 
lution of the outburst (see also iKolotilov fc Kenvon1ll997f) . 
We discuss below in section|l]our wind-envelope interaction 
model to explain this recent luminosity drop. 

3.2 The colour evolution of V1515 Cyg 

As described above, the long rise time of light curve of 
V1515 Cyg suggests that in this case the thermal instabil- 
ity is_jirjt_t£iggered, but proceeds inside-out, as described 
bv lBell fc Linl (fT994Tl . We have constructed time dependent 
models also in this case and compared them to the colour 
evolution observed in this system. In this case, we have 
adopted M* = IMq and R m i n = 3Rq. We have also as- 
sumed Menv = 2 x 10 -6 MQ/yr. With these, parameters, we 
obtained a slow-rise outburst, with a peak luminosity of 250 
Lq and a peak mass accretion rate M pca k ~ 1 x 10 _4 MQ/yr. 
Fig.|H|shows the colour evolution in the V — (B — V) plane 
in this case. V1515 Cyg is assumed to be at a distance of 
1000 pc, with a face-on disc. The models have been red- 
dened assuming Ay = 3.1 and R = 3.3. The assumed value 
of Ay is slightly smaller than the one used for V10 57 Cyg, 
consistent with previous estimates jBell et alJll99oTl . As be- 
fore, the solid line shows the result from our time dependent 
model, while the dashed line show the results from a series 
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of steady state models. The agreement is quite good, also 
given the little evolution observed in V1515 Cyg. 



4 INTERACTION BETWEEN DISC WIND 
AND INFALLING ENVELOPE 

4.1 Evidence for wind and envelope components 
in FU Orionis systems 

As mentioned above, it has long been recognised that ad- 
ditional luminosity components are required in order to ex- 
plain a range of spectral diagnostics. Here we focus on an 
infall ing dusty envelope (invoked by iKenvon fc HartmanrJ 
Il99ll in order to explain the magnitude of the 10/xm ex- 
cess, particularly in V1057 Cyg and V1515 Cyg) and the 
strong winds in FU Orionis systems (as evidenced by P 
Cygni profiles and 'shell' features in optical spectra) . Before 
proceeding to consider the dynamical interaction between 
these components, we briefly describe the observational ev- 
ide nce for (and the inferred properties of) each. 

IKH91I introduced the notion of an infalling envelope in 
order to explain the magnitude of the far infrared flux (i.e. 
at A > 10/xm) in FU Orionis systems, which is underpre- 
dicted by steady state, flat accretion disc models. Although 
several other suggestions have been made regarding the ori- 
gin of the fa r infrared emission in FU Orionis objects (see , 
for example. lLodato fc Bertinll200ll : lAbraham et alJl2004> . 
the similar fading rates of the optical emission and that in 
the range 10 — 25/im suggests that the latter is derived from 
reprocessing of optical light (although the near constancy of 
emission longward of ~ 25/im suggests an additional compo- 
nent must supply the longest w avelength emission, see dis- 
cussion in lAbraham et alll2004h . The dusty envelope how- 
ever remains the most likely site of emission in the 10 — 25pim 
range, since alternative scenarios (involving reprocessing in 
a flared accretion disc) are hard to construct, given the ex- 
pected curvature of the disc surface in realistic disc models 
ilBell et al.lfl997l) . 

The spectral modeling of IKH9 3 implied that the en- 
velope must subtend a reasonably large solid angle at the 
source (covering factor ~ 0.5) but that the observer views 
the source at rather low inclination and thus sees a low opti- 
cal extinction to the source. This therefore implies that the 
envelope is somewhat flattened (but not disc like) consistent 
with the ex pected morphology for a rotating infalling enve- 
lope. EH<ni assumed a density profile appropriate to free fall 
onto a point mass 



M c , 



47rr-o /2 (2GM*) 1 /2 V r > 



(2) 



and adjusted the normalisation (and hence infall rate in the 
envelope) in order that the flux reprocessed in the optically 
thin dust envelope matched the spectral energy distribu- 
tion longward of lO/xm. An attractive feature of this model- 
ing was that the inferred infall rate (4 10 _6 Mq yr" 1 ) is of 
the correct magnitude to trigger thermal instability in the 
inner disc (see Section l3.1|l . However, since the disc mod- 
els can themselves provide the observed flux shortward of 
lO^im (and are c onsistent with the tem poral behaviour of 
these wavebands; lAbraham et alJl2004h . the envelope must 
not make a significant contribution shortward of ~ 10/^m. 



This requirement constrains the inner edge of the envelope 
to lie at around 10 au from the central source. Since this 
radius is considerably outside the dust sublimation radius 
(which, for a central source luminosity of a few hundreds 
1/0, lies at a distance of « 1 au) this begs the question of 
what prevents the envelope extending to smaller radii. 

The presence of the wind is unambiguously indicated by 
P Cygni profile s in Ca II A8542, Na I and the Balmer lines o f 
hydrogen ijCroswell et aljfl987l : lHartmann fc Calvetl [l995). 
indicating outflow velocities of up to ~ 400 km s _1 . Like- 
wise, lower excitation lines of neutral metals and TiO bands 
exhibit 'shell features', i.e. absorption features displaced by 
about 50 km s _1 to blueward, which are inter preted as aris- 
ing from localised condensa tions in the wind llHerbig et alJ 
2003.: lHartmann et al.l2004h . Detailed modeling of these fea- 
tures yield wind mass loss rates which correlate positively 
with the inferred accretion rate onto the star, as is observe d 
in other classes of young stars l|Hartmann fc Calvetl ll995V 

4.2 Wind-envelope interaction 

The evolution of a wind expanding into a quasi-spherical 
densi ty distribution is initially approximately adiabatic 
fe.g.. IWeaver et ai1ll977l) . i.e. the combined thermal energy 
of the wind blown cavity and the kinetic energy of the swept 
up shell of shocked medium together increase at a rate set 
by the rate of kinetic energy input in the wind (~ M w wJ, 
where u w is the terminal velocity of the wind). In a flat- 
tened density distribution, however, the wind 'breaks out' 
in the polar direction. We simplify the problem by consider- 
ing a wind that expands into a free-falling envelope, whose 
density (in spherical polar coordinates) is given by equa- 
tion ^ for 8 > 9 C , where C measures the angular extent of 
the circumpolar cavity (see Fig. for a schematic view of 
the assumed geometry). In this case, the envelope absorbs 
a fraction / = cos C of the momentum of the wind. If the 
rate of momentum input in the wind is p = /M w « w , then 
the total momentum provided to the flow at time t is given 
by 



p(*)= fp{t')dt'. 
Jo 



(3) 



The wind drives a strong shock into the envelope, of ra- 
dius R s {t) and the gas swept up by the shock is mainly con- 
centrated in a thin shell at ~ R s . For the density distribution 
described in Eq. (|5J and assuming that the envelope is ini- 
tially in free fall (i.e with inward velocity ur = (2GM*/ R) 1 ^ 2 
for central object of mass M*), the mass contained within 
R s is given by: 

1/2 



M(R S ) = 



( 2Rj 
{ GM* 



(4) 



In the limit that the gravitational acceleration of the 
flow is small compared with that resulting from interaction 
with the wind, the evolution of a momentum conserving thin 
shell is given by the condition: 

p{t) = M(R S )R S . (5) 

In the simple case that p(t) is a simple power law 
p(t)=p (t/t ) a , (6) 
we can readily obtain power law solutions for R s , i.e. 
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Figure 9. Schematic view of the system geometry assumed here. 



R s 



Ro(t/t ) 



where 



Ro = 



V2 3 / 2 (l + a)Af cl w 



2/r, 



(7) 



(8) 



and b = 2/5(o + l). 

Therefore any declining rate of momentum input im- 
plies an evolution intermediate between the constant mo- 
mentum (R B oc t 0A ) and the constant p (R B oc i ' 8 ) case. 

The neglect of the thermal and ram pressure of the 
material over-run by the shock implies that this solution 
is valid only if _R S is greater than both the sound speed, 
c s and the local free fall velocity of the unshocked enve- 
lope, wr. Since the envelope is initially flowing in super- 
sonically, the latter criterion is more restrictive. The sim- 
ilarity solution derived above implies R B oc /j2a-3/2(a+i)^ 
whereas vr oc R^ 1 ^ 2 . Thus R s /vn declines with increasing 
R s provided a < 2/3. Therefore whereas for a steady wind 
(a = 1), gravity becomes of decreasing importance as the 
shock propagates outwards, in the case of a one-off appli- 
cation of momentum (a = 0), the shell will be increasingly 
subject to gravitational retardation and may be expected to 
re-collapse. We shall be applying a rate of momentum input 
from the wind that tracks the luminosity evolution of the 
system. In the case of FU Orionis, we adopt a = 1, appro- 
priate to the nearly constant system luminosity, whereas in 
V1057 Cygni, the exponential decline in system luminosity 
means that at late times, we expect the behaviour to ap- 
proach that of a one-off impulse. We therefore expect the 
dust shell to expand monotonically in the case of FU Orio- 
nis and for it to re-collapse in the case of V1057 Cygni. In 
the latter case we may estimate the impulse applied by the 
wind (~ M w i; w At), where At is a measure of the effective 
timescale over which the bulk of the momentum is applied 
(~ a year for this rapidly fading system). Equating this to 



p(t), we can therefore obtain a rough estimate of the radius 
to which the shell should be swept out (equations @-(|SJl). 
We get: 



Rou 



3/ M, 



(9) 



If we adopt the parameters we use for the simulations of 
V1057 Cygni below (i.e. M„ = 6.3 x l(T 6 M /yr, M cnv = 
10 Mq/vt, / = 0.5 and u w = 300 km/sec), we estimate 
that the shell should be swept out to radii ~ 30 au before 
re-collapsing. 



4.3 Simulations 

We explore the interaction between the time dependent mo- 
mentum input from the wind and the infalling dusty enve- 
lope using a simple one dimensional L agrangian code which 
empl oys a standard artificial viscosity jRichtmver fc Morton! 
1967) in order to model the shock at the interface between 
the wind and envelope. For modeling the evolution of V1057 
Cygni, we adopt the parameters listed at the end of 4.2, and 
assume a stellar mass of M* = 0.5Mq. In the absence of 
a proper calculation of the gas temperature in such an en- 
velope, we simply model the flow as isothermal gas with 
temperature 600K. A temperature of hundreds of Kelvin is 
a reasonable order of magnitude etsimate based on the ex- 
pected temperature of large grains at distances of 1 — 10 au 
from a source with luminosity of hun dreds of Lp, and issim- 
ilar to the temperature employed bv lMuzerolle et ail <l2005f) 
in modeling the SED of the recently erupted FU Orionis 
system V1647 Ori. We note that the temperature does not 
appear in the equations above for the evolution of a thin 
shell, since the only requirement here is that the shock is 
strong, a condition that is readily met (with this choice of 
temperature) provided the shock remains within « 100 au. 
We will show below that doubling the temperature changes 
the outermost propagation radius and associated timescale 
by less than 30%, so our conclusions about the applicability 
of this model to FU Orionis systems are not unduly sensitive 
to our crude treatment of the thermal structure of the flow. 

The impulse is initially applied to fluid elements out- 
ward of radius -Rin, smoothed over a half Gaussian with 
standard deviation a = 0.1-Rm- Since we expect the enve- 
lope to intercept the disc at small radius, R[ n is a measure 
of the radius at which the envelope starts to subtend a sig- 
nificant solid angle as seen from the source, and therefore 
corresponds to the region where the initial wind-envelope 
interaction occurs. In what follows, we adopt, somewhat ar- 
bitrarily, Ri n — 2 au. We note that the actual value of Ri n 
is not critical in determining the late time evolution of the 
system, unless it is so small that the ram pressure of the 
envelope exceeds the initial momentum input rate from the 
wind. In practice, for typical parameters for the winds in 
these systems, a quasi-spherical envelope would have to ex- 
tend to radii < 0.1 au in order to be able to crush the wind 
in this way. 

Fluid elements inwards of _Ri n are rapidly accreted, so 
that the density distribution develops an inner cavity that 
propagates outwards behind the shock. Provided that fluid 
elements initially at radii > Ri n remain in this region, the 
wind momentum continues to be applied to them according 
to the initial (Gaussian) distribution of relative weightings. 
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Figure 10. Snapshots of density evolution profile in the envelope. The left hand panel depicts the outward propagation of the swept up 
shell (at times 4, 8 and 12 years after the onset of outburst for successively larger outer shell radii), while the right hand panel shows the 
shell re-collapse at 16, 20 and 24 years. 



If an element initially at radius > i?i n subsequently falls 
within Ri n , it is no longer subject to momentum input from 
the wind, and the set of relative weightings are passed up- 
stream by one fluid element. In this way, momentum input 
is always distributed in the same way over the fluid elements 
that are currently at radius > Ri n . 

For the overall normalisation of the wind momen- 
tum input rate we assume that this varies in propor- 
tion to the total luminosity of the source, motivated by 
the general correspondenc e between the two discussed by 
lHartmann fc Calved lll995l) . Note that although we here as- 
sume that this also roughly applies to the evolution of in- 
dividual sources, this has not been demonstrated to be the 
case: although V1057 Cygni, with its dramatic fading over 
its first decade in outburst, would have been a good source 
with which to test this hypothesis, detailed wind model- 
ing was only performed in the late 1980s, when the sys- 
tem had already faded by a large factor. Here we model 
the wind of FU Orionis as a constant mass input rate of 
6.3 x lO~ 5 M /yr and that of V1057 Cy gni as having a peak 
mass input rate of 6.3 x 10 -6 MQ/yr, in both cases assuming 
a wind terminal velocity of 300 km/s. We parameterise the 
decline of the wind in V1057 Cygni as an exponential de- 
cline (timescale 6.4 years) for the first 4.5 years followed by a 
slower exponential decline (timescale 15.5 years), thereafter. 

In order for the results to be independent of resolution, 
we find that it is necessary for no more than 20% of the 
total wind momentum to be applied to a single grid cell. In 
practice, this implies AR/R < 0.03 at the inner edge; in our 
calculations, the shell is driven out to radii 10 — 100 au and 
we need to model the flow out to ~ 50 au in order that the 
front propagation is unaffected by the outer boundary of the 
flow. We ran convergence test on our simulations, employing 
up to 24000 equal mass elements. 



4.4 Results 

Figure [TTJ1 illustrates snapshots of the density profile of the 
envelope for the V1057 Cygni models whose parameters are 
listed above. As expected, the flow develops an internal cav- 
ity bounded by a thin layer of shocked material, whose den- 
sity exceeds that of the unshocked infalling envelope by the 
large factor (~ the square of the Mach number) expected 
for strong isothermal shocks. As expected, given the steep 
decline in the momentum input from the wind, the shock 
weakens as it propagates outwards and the density contrast 
at the shock consequently decreases. Gravitational decelera- 
tion causes material in the shell to - more or less coherently 
- reverse direction about a decade after outburst, when it 
has reached a radius of w 10 au. The shell-like structure 
is initially preserved during the fallback stage (see the den- 
sity peak at ~ 6 au in the right hand panel of Figure 1101 
at 16 years after outburst); the shell is initially infalling at 
less than the free fall velocity and is therefore overtaken by 
faster material flowing in from the envelope. As the shell 
picks up speed, due to gravitational acceleration, the veloc- 
ity differential with respect to the inflowing envelope suc- 
cessively decreases. Once this differential becomes subsonic, 
the shock disappears and the shell-like structure is eroded, 
as material at its inner edge is acclerated ahead of the rest 
of the flow (the remnant shell at ~ 2 au is visible in the 
snapshot at 20 years in the right hand panel of Figure HOB . 
By the stage that the flow has fallen back to small radii, the 
density profile has relaxed to a smooth p oc r~ 5 for free 
fall onto a point mass. 

We also ran a simulation identical to that above, except 
with a temperature of twice the above value (i.e. 1200K com- 
pared with 600K). As one would expect, the evolution was 
qualitatively very similar during the outward propagation of 
the shell (Figure mi . but differences appeared close to turn- 
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Figure 11. Snapshots of density evolution profile for the V1057 Cygni simulation with T = 1200K. The left hand panel depicts the 
outward propagation of the swept up shell (at times 4, 8, 12 and 16 years after the onset of the outburst for successively larger outer 
shell radii), while the right hand panel shows the shell re-collapse at 20, 24 and 28 years. 



round, where the shell motion is temporarily subsonic. As 
expected, the shell is considerably thicker in the case of the 
hotter simulation, and by the time the velocity of the entire 
shell is reversed, the shell thickness is already 50% of its ra- 
dius (see right hand snapshot in the right panel of Figure 
1111 . The re-collapse never proceeds by way of a thin shell (as 
in the colder simulation); instead, the flow passes through a 
shock that remains more or less stationary in the frame of 
the star on a timescale of a decade or so: the inward flowing 
material is decelerated in the shock and then gravitationally 
accelerated, finally free-falling onto the star. In contrast to 
the colder simulation, a pronounced shock structure remains 
at round 10 au, even « 30 years after the onset of outburst. 

In Figure rHn by contrast, we show how the shock radius 
evolves for our (steady wind) model for FU Orionis. Here, 
(as shown by the dashed line in Figurc tT^l . the radius of the 
shock evolves according to the R s (xl 0,8 similarity solution 
derived in equation J7J above. As we noted above, R„/vr 
increases with time, so that the role of gravity becomes less 
significant as the front propagates outwards. At the current 
epoch, we would expect the shock to have propagated to 
~ 1000 au; if FU Orionis remained in outburst indefinitely, 
then the shock would reach the outer boundary of the star's 
putative natal core (at ~ 10 4 au) after about 10 3 years, and 
would still be strongly supersonic at this point. 



log t(years) 

Figure 12. The evolution of the radius of the swept up shell for 
the steady wind simulation for FU Orionis. For comparison, the 
dashed line depicts the R a oc t°- s similarity solution derived in 
equation 171 . 



4.5 Discussion 

We have modeled the interaction between the wind and in- 
falling envelope of FU Orionis systems using parameters for 
these components that have previously been derived in the 
literature. This modeling has demonstrated that there are 
several consequences of this interaction that can explain a 



number of characteristics of FU Orionis systems. In partic- 
ular, we are tempted to make the following connections: 

i) that the inner radius of the dust envelope (invoked in 
order to generate the 10 — 15^im flux of certain FU Orionis 
systems via reprocessing of the radiation from the central 
source) is set by the dynamical interaction between the wind 
and the inflowing envelope 
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and ii) that the erratic photometric variability (ob- 
served in V1057 Cygni post 1996 and also in V1515 Cygni) is 
associated with the fall back of dusty material to small radii 
and the consequent passage of dust condensations across the 
line of sight to the inner accretion disc. 

Below we examine each of these hypotheses in turn: 

4-5.1 Effect on the mid infrared SED 

Figure 1101 shows that in our model for V1057 Cygni, the 
wind sweeps out a cavity that attains a size of 10 au about 
five years after the onset of the outburst, and which remains 
in the range 10 — 20 au for the next fifteen years or so. A 
cavity of this dimension is precisely what was invoked by 
IKH91I in their modeling of the mid infrared excess of V1057 
Cygni, based on observations made 10 — 20 years after the 
outburst. Our model is therefore perfectly compatible with 
these results. Thereafter, in our model, material falls back to 
small radii, although, depending on the temperature of the 
gas, the density discontinuity at the furthest extent of the 
windblown bubble (at about 10 au) may persist for decades 
after that. In the absence of two dimensional modeling, both 
of the hydrodynamic interaction and of transfer of radia- 
tion through the resulting structure, it is unclear to what 
extent we would expect this relic structure to leave an im- 
print on the spectral energy distribution. We therefore do 
not attempt to predict how the SED would evolve during 
this fallback phase, pending more detailed calculations. In 
our models for FU Orionis, by contrast, the wind has swept 
the envelope to radii « 10 3 au and we would not expect to 
observe a spectral signature of a reprocessing cavity in the 
mid infrared. Again, this is compatible with observations of 
FU Orionis, where the evidence for extra spectral compo- 
nents in the mid infrared (over and above what is supplied 
by a f iat acc retion disc) is considerably weaker than in V1057 
Cyg JKHU). 

Although we have not modeled V1515 Cyg in any de- 
tail, we note that the weaker wind of this system (Hartmann 
and Calvet 1995) leads one to expect that the envelope is 
not expelled to very large radii. In this respect one would 
expect the system to more resemble V1057 Cyg than FU 
Orionis. In fact, V1515 Cyg is indeed similar to V1057 Cyg 
to the extent that it exhibits a significant mid infrared ex- 
cess. Finally, we note that in the recently erupted rapid rise 
system V1647 Ori (McNeil's nebula) , modeling of the SED 
derived from Spitzer and ground based data again required 
an infalling envelope with infall rate of ~ 1O -6 M0 yr _1 
jMuzerolle et al.ll2005h . In this modeling, the inner edge of 
this envelope was placed at 1 au; this compact inner edge 
is expected in this model a the early stage of the outburst 
at which this system was modeled; depending on the sub- 
sequent luminosity evolution during the outburst, we might 
expect this radius to propagate outwards over the next few 
years. 

4-5.2 Effect on photometric variability 

Figure 1101 shows that after around 20 years the dusty en- 
velope in our model for V1057 Cygni has collapsed back 
to small radii and it is obviously tempting to associate 
this event with the sudden dimming of the system in 1996 



an d its subsequent erratic va riability. (Note that, follow- 
ing |^^i^TO^^^^k^3Ei3, we can readily demonstrate 
that, given the densities encountered in these simulations, 
the dust should not be destroyed in the shock and should 
instead melt only when it falls back within the dust sub- 
limation radius at less than an A.U.). As has been noted 
(KH91; see also Figure [7] Section 3) the colour variations 
of these fluctuations are consistent with screening of the 
source with a variable dust screen following the standard 
interstellar extinction law (Mathis 1990). It is however nec- 
essary to postulate that in the realistic flow morphology for 
a recollapsing, rotating envelope, most of the envelope joins 
the disc at radii > a fews tenths of an au, so that only a 
small fraction of its mass can intercept the line of sight to 
the innermost disc which produces the bulk of the optical 
emission. If this were not the case (i.e. if - as in our simpli- 
fied one dimensional model - the flow just collapsed back, 
re- filling a wedge of constant opening angle), the extinction 
through the envelope would be immense (Av > 1000) and 
under these circumstances the dimming of the source would 
be much more dramatic than that observed (~ 1 magni- 
tude). In fact, it is a quite realistic expectation that only 
a small fraction of the envelope material should be on such 
low angular momentum trajectories that it falls in to such 
small radii, but without a two dimensional model we cannot 
quantify this further. 

Such a pictu r e in s ome ways resembles that proposed 
bv lKenvon et all lll99ll) for the similar photometric varia- 
tions of V1515 Cygni a nd also for the more recent behaviour 
of V1057 Cygni itself jKolotilov fc Kenvodll997h . In both 
cases, the variations derive from variable extinction local to 
the source. However, Kenyon et al and Kolotilov and Kenyon 
postulate that this dust is formed afresh in the wind, in a 
manner analogous to that observed in the winds of classical 
novae iGehrzll98al . rather than belonging to the re-collapse 
of a dusty accretion flow. If the dust is co-moving with the 
wind, then the timescale for the observed variability (of or- 
der a year, see section l2~2l requires that, in t heir model, the 
dust be situated at large radii (~ 100 au, see lKenvon et all 
1991), whereas in our model the dust is roughly in a state of 
free fall and therefore must be located at around au distance 
scales. 

There are several points in favour of this interpreta- 
tion. Firstly, the fallback of the envelope provides an ex- 
planation of the change in the photometric properties of 
V1057 Cygni after 1995. It is also consistent with the pho- 
tometric variability patterns of other FU Orionis systems. 
V1515 Cygni has demonstrated similar photometric vari- 
ability patterns to those recently observed in V1057 Cygni 
throughout its outburst. The wind in V1515 Cyg is rela- 
tively weak, comparable to that measured in V1057 Cygni 
after its had undergone significant fadin g from the peak of 
the outburst IIHartmann fc C alvet 1995). We therefore pos- 
tulate that whereas V1057 Cygni was able to temporarily 
clear its circumstellar environment in response to the strong 
burst of wind activity at the onset of its outburst, the wind 
in V1515 Cygni has never succeeded in completely clearing 
the line of sight to the observer of dusty material and hence 
the system has been subject to erratic variations throughout 
the outburst. In FU Orionis, by contrast, the stronger, more 
sustained, wind has succeeded in clearing a large cavity and 
there is no trace in its light curve of the sort of variability 
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exhibited in the other two systems. If, on the other hand, 
the variations result from dust formed in the wind, there is 
no particularly obvious reason why the wind in FU Orionis 
has apparently never undergone dust condensation events, 
that in V1515 Cygni has been subject to dust condensation 
events throughout and the wind in V1057 Cygni has made 
an abrupt transition between FU Orionis like and V1515 
Cygni like behaviour. 

The other advantage of our interpretation is that it does 
not predict any straightforward relationship between pho- 
tometric variations and the appearance of prominent shell 
structures in the optical and infrared spectra. These shell 
structures correspond t o blue shifted absorption features at 
aroun d 50 — 100 km/s llHerbig et al.ll20o3 : lHartmann et alJ 
l2004f) and are interpreted as arising from shell-like regions 
of enhanced density in the outflowing wind. There has been 
some discussion of the hypothesis that the shell could be 
the formation site of the dust associa t ed with the photo- 
metric variations teenvon et all 11991k iHerbig et alTl2003l 
lHartmann et alJl2004f) . an idea given additional credence by 
the fact that the shell features in V1057 Cygni became very 
strong when the system faded dramatically at the end of 
its plateau phase. (We note in passing that the dust forma- 
tion would have to be very inefficient in this case to account 
for the amplitude of such variatio ns: the inferred colum n 
density in the shell - ~ 10 23 cm" 2 : lHartmann et"ai]|2004l - 
would correspond to 100 magnitudes of optical extinction 
if the dust formed with a standard dust-gas ratio). More 
importantly, however, the shell feature would be expected 
to be at its strongest at minimum light of the system (in 
1999) whereas, in fact, the shell features had essentially dis- 
appeared at this date, to re-appear a couple of years later 
as the system recov ered in photometric brightness. This led 
IHerbig et alJ ((2003) to suggest that the photometric varia- 
tions instead result from 'inhomogeneities ... passing across 
the line of sight', a situation that would result from the sort 
of clumpy fall back that we envisage here. In our model, al- 
though the enhanced shell features in 1996 may or may not 
have anything to do with the fall back of the envelope, the 
dust is not generated in, nor is necessarily cospatial with, 
the shells and one would not expect to be able to trace the 
effect of individual shell ejection events in the light curve of 
the object. 

The largest uncertainties in what we have sketched 
above concern how the returning flow would behave in two 
dimensions and what is the degree of clumping that would 
be expected in reality. Two dimensional modeling would 
remedy the first shortcoming: one might readily anticipate 
that Ray leigh- Taylor instabilities in the interface between 
the infalling envelope and wind would produce a wealth of 
time dependent behaviour that could be checked against the 
light curves of V1515 Cygni and V1057 Cygni. However, one 
would also expect that much of the fine grained structure in 
the light curves of these objects would result from small 
scale clumping that would not be resolved in such simu- 
lations. Certainly, HST imaging of V1057 Cygni reveals a 
wealth of features (loops and arcs) in the dust distribution 
on a larger scale (hundreds of au), and it does not seem un- 
reasonable to suggest that the material falling in to small 
radii would be far from homogeneous. In the absence of 
a more quantitative argument, the best observational test 
of our hypothesis remains the continued spectroscopic and 



photometric monitoring of V1057 Cygni and V1515 Cygni 
in order that our prediction (of no correlation between in- 
dividual shell ejections and the photometric variability) can 
be properly assessed. 



5 CONCLUSIONS 

Understanding the FU Orionis phenomenon is an important 
issue in the general context of star formation, especially if, as 
is commonly thought, most protostars undergo an FU Ori- 
onis phase during their early evolution. In particular, given 
their complexity and the different components involved in 
the phenomenon (protostar, accretion disc, outflows, enve- 
lope, possible small mass companions) their represent a valu- 
able tool to study the interaction of the protostar with its 
environment. 

Here, we present new UBVR photometric data of the 
three best studied FU Orionis objects (FU Ori, V1057 Cyg 
and V1515 Cyg), taken at Maidanak Observatory between 
1981 and 2003. We have analysed our photometric data for 
variability on long and short timescales. We have therefore 
been able to monitor the recent evolution of these systems. 
Whereas V1515 Cyg and FU Ori did not show any signif- 
icant change in their long-term evolution, V1057 Cyg has 
undergone an abrupt dro p in luminosity in 1995 (see also 
iKolotilov fe Kenvonlll99'7l) . Apart from monitoring the sec- 
ular evolution of these systems, we have also checked for the 
presence of variability on shorter timescales and have de- 
tected non periodic, flickering, behaviour on a timescale of 
a few days both in V1515 C yg and in V1057 Cy g, analogous 
to that observed in FU Ori llCenvon et aLE oOP). Both these 
objects also displayed periodic variations (period ~ 14 days) 
during one observing season in each case. 

We have then interpreted the secular evolution of these 
systems within the framework of two complementary the- 
oretical models, both of which involving the interaction of 
different components of the system. 

(i) The overall outburst and subsequent luminosity de- 
cline (especially for V1057 Cyg, that displays the fastest 
evolution) has been described in terms of a standard disc 
instability model (the instability being possibly triggered by 
the interaction between the circumstellar disc and a massive 
planet or a small mass companion). For the first time, we 
have compared photometric data at different epochs with 
time-dependent disc models, rath er than with a series of 
steady state models (as done by KH91), finding a better 
agreement, with respect to the latter, of the colour evolu- 
tion of the systems, particularly in the V — (B — V) colour- 
magnitude diagram. In particular, we have shown that we 
can reproduce the large change in system colour between 
1970 and 1995 which cannot be produced by a sequence of 
steady state disc models (see Figure 0. 

(ii) On the other hand, we have interpreted the strong 
luminosity drops observed early during the outburst in 
V1515 Cyg and more recently also in V1057 Cyg as oc- 
cultation effects due to intervening dust along our line 
of sigh t. A similar interpretation has already been pro- 
posed jKenvon et alJll99lt Ikolotilov fc Kenvonlll997l) . but 
whereas these authors assume that the dust is formed within 
an expanding shell in the wind, here we consider the occul- 
tation events as the result of the interaction between the 
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wind and the surrounding infalling envelope. We have con- 
structed a simple quasi-spherical model of the evolution of 
a dusty envelope subject to the momentum input from the 
outflowing wind. For sufficiently strong winds, the envelope 
is blown out to large radii (w fO au, consistent with the 
required location of the reprocessing dust, in order to pro- 
duce t he mid-infrared excess observed in some systems, see 
IKH91I) . However, at later stages during the outburst, as the 
momentum carried by the wind decreases (as is expected if 
the wind strength scales with the luminosity of the disc), the 
dusty envelope falls back and occults our line of sight. In this 
way, our model naturally explains the different photometric 
variability of the three FU Orionis objects studied here in 
terms of the different strength and time dependence of the 
disc wind. Objects with strongest winds (FU Ori) are able 
to clear away the envelope and do not show any significant 
luminosity drop, those with weaker winds (Vf 5f 5 Cyg) show 
erratic photometric variability throughout , whereas Vf 057 
Cyg, which has faded substantially in the last 20 years, has 
started to show variability only after the (disc) luminos- 
ity has decreased substantially. The main limitation of our 
model is its fD nature, which makes it difficult to quantify 
in greater detail the amount of occulting material that will 
intercept our line of sight after the envelope has fallen back 
to small radii. We leave this issue to further investigations. 
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